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Abstract 


Improving  the  performance  of  most  engineering  systems  requires  the  ability  to  model  the  system's  behavior 
with  improved  accuracy,  'llic  evolution  of  the  mechanical  arm  from  tclcopcrator  and  crane  to  present  day 
industrial  and  space  robots  and  large  space  manipulators  is  no  exception.  Initial  simple  kinematic  and 
dynamic  models  are  no  longer  adequate  to  improve  performance  in  the  most  critical  applications.  Both  the 
mechanical  system  and  control  system  require  improved  models  for  design  simulation.  Proposed  new  control 
algorithms  require  dynamic  models  for  control  calculation.  Planning  and  programming  activities  as  well  as 
man-in-thc*loop  simulation  also  require  accurate  models  of  the  arms. 

Accuracy  is  usually  acquired  at  some  cost  The  application  of  mechanical  arms  to  economically  sensitive 
endeavors  in  industry  and  space  also  gives  incentive  to  improve  the  efficiency  of  the  formulation  and 
simulation  of  dynamic  models.  Control  algorithms  and  man-in-the-loop  simulation  require  real  time 
calculation  of  dynamic  behavior.  Formulation  of  the  dynamics  in  an  easy  to  understand  conceptual  approach 
is  also  important  if  maximum  use  of  the  results  is  to  be  obtained. 

The  nonlinear  equations  of  motion  for  flexible  manipulator  arms  consisting  of  rotary  joints  connecting  two 
flexible  links  are  developed.  Kinematics  of  both  the  rotary  joint  motion  and  the  link  deformation  are 
described  by  4x4  transformation  matrices.  The  link  deflection  is  assumed  small  so  that  the  link 
transformation  can  be  composed  of  summations  of  assumed  link  shapes.  The  resulting  equations  are 
presented  as  scalar  and  4x4  matrix  operations  ready  for  programming^The  efficiency  of  this  formulation  is 
compared  to  rigid  link  cases  reported  in  the  literature. 

Keywords:  Robots;  Distributed  parameter  systems;  Models;  Manipulation;  Vibration  control;  Flexible 
mechanisms;  Mechanical  arms. 
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1.  Sketch  of  Prior  Work 

Much  work  has  been  done  to  formulate  the  dynamic  equations  of  motion  for  mechanical  arms  with 
rigid  links.  Work  on  the  "inverse  dynamic  formulation"  used  in  control  can  be  found  in  references  [22],  [27], 
[29],  [2]  and  in  their  bibliographies.  References  [301.  [20],  [33],  [32],  [12],  and  their  bibliographies  represent 
work  on  the  dynamic  formulation  for  simulating  rigid  link  arms.  The  efficiency  of  these  formulations  and 
alternatives  to  their  real  time  calculation  is  discussed  in  [26],  [1]  and  the  works  referenced  therein. 

'Che  limitation  of  these  works  is  that  rigid  links  arc  assumed.  With  this  assumption  the  techniques 
become  at  some  point  self  defeating,  if  their  purpose  is  to  improve  performance.  Maintaining  rigidity  of  the 
links  inhibits  improved  performance  but  is  necessary  if  the  rigid  link  assumption  is  to  be  accurate. 

Consideration  of  flexibility  and  control  of  the  links  in  arm-type  devices  was  reported  in  1972  by 
Mirro  [24].  Ihis  early  work  considered  both  the  modeling  and  control  of  a  single  link  device.  Book  [7] 
considered  the  linear  dynamics  of  spatial  flexible  arms  represented  as  lumped  mass  and  spring  components 
via  4x4  transformation  matrices.  This  was  refined  and  later  reported  in  [9|.  Rook  and  Whitney  [3],  [4]  later 
considered  linear  distributed  dynamics  of  planar  arms  via  transfer  matrices  and  the  limitations  flexibility 
imposed  on  control  system  performance  [8].  Maizza  and  Whitney  [23],  [4]  used  a  planar  nonlinear  model  with 
modal  representation  of  the  flexibility  and  considered  modal  control  as  a  technique  for  overcoming  the 
limitations  of  the  flexibility.  Whitney,  Book,  and  I  .ynch  [34],  [4]  considered  the  design  implications  of 
flexibility.  Distributed  frequency  domain  analysis  of  nonplanar  arms  using  transfer  matrix  techniques  [S],  [6] 
has  been  used  by  Book,  ct.al  to  verify  the  accuracy  of  truncated  modal  models  of  the  nonlinear  spatial 
dynamics  of  flexible  manipulators  (the  Remote  Manipulator  of  the  Space  Shuttle).  The  nonlinear  modal 
model  appearing  here  was  first  presented  by  the  author  in  1982  [10].  A  more  classical  approach  to 
manipulator  dynamics,  both  rigid  [18]  and  flexible  [19],  has  been  undertaken  by  Huston  and  his  coworkers. 

Ihc  work  in  flexible  spacecraft  has  spawned  a  line  of  research  pertaining  to  the  interaction  of  articulated 
structures.  This  work  has  great  relevance  to  the  manipulator  modeling  problem.  Kntrics  into  this  literature 
arc  provided  by  the  works  of  1  .ikins  [21]  and  Hughes  [25].  This  activity  produced  a  spatial,  nonlinear,  flexible 
manipulator  model  reported  by  Ho  eual.  [14]  and  corresponding  computer  code  for  simulation.  The 
simulation  required  great  amounts  of  computer  time  and  was  unsuitable  for  even  offline  simulation.  Further 
work  for  the  purposes  of  simulating  the  Space  Shuttle  Remote  Manipulator  was  performed  by  Hughes.  His 
linearized  model  is  reported  in  [16]  and  a  more  general  model  is  reported  in  [17].  The  Hughes  model  ignores 
the  interaction  between  structural  deformation  and  angular  rate  as  might  be  appropriate  for  the  Space  Shuttle 
arm.  This  work  and  associated  work  at  SPAR  Aerospace.  Ltd.  and  the  Charles  Stark  Draper  laboratory.  Inc. 
probably  represent  the  most  intensive  work  on  the  modeling,  simulation  and  control  of  flexible  arms. 
Unfortunately,  little  of  this  work  has  been  reported  in  the  open  literature.  Recent  examination  of 
experimental  results  from  the  operation  of  the  Shuttle  arm  in  space  has  confirmed  the  validity  of  these 
models.  More  recently.  Singh  and  Likins  [28]  have  reported  an  efficient  flexible  arm  simulation  program. 

Yet  another  branch  of  research  that  has  found  its  way  to  the  flexible  manipulator  dynamics  problem  is 
the  study  of  flexible  mechanisms.  Dubowsky  and  Gardner  [13]  and  Winfrey  [35]  provide  the  reader  with  a 
bibliography  on  this  work.  Sunada  and  Dubowsky  [31]  have  developed  modeling  techniques  applicable  to 
both  spatial  dosed  loop  mechanisms  and  open  loop  chains  such  as  manipulator  arms,  litis  work  assumes  a 
known  nominal  motion  over  time  about  which  die  flexible  arm  equations  arc  linearized.  This  fafls  short  of  a 
true  simulation  of  the  flexible,  nonlinear  equations,  but  is  an  interesting  compromise  for  the  sake  of 
computational  speed.  This  technique  is  oriented  toward  finite  element  analysis  to  obtain  modal 
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characteristics  of  the  links  which  arc  then  combined  using  a  time  varying  compatibility  matrix.  It  uses  4x4 
matrices  to  represent  the  nominal  kinematics  and  derivation  of  the  compatibility  matrix. 

1.1.  Perspective  on  This  Work 

'Ibis  report  stresses  an  efficient,  complete,  and  conceptually  straightforward  modeling  approach  using 
the  4x4  transformation  matrices  that  arc  familiar  to  workers  in  the  field  of  robotics.  It  is  unique  in  several 
respects:  It  uses  4x4  matrices  to  represent  both  the  joint  and  deflection  motion.  Ihc  deflection 
transformation  is  represented  in  terms  of  a  summation  of  modal  shapes.  I  hc  computations  resulting  from  the 
lagrangian  formulation  of  the  dynamics  arc  reduced  to  recursive  form  similar  to  that  which  has  proven  so 
efficient  in  the  rigid  link  ease.  The  equations  arc  free  from  assumptions  of  a  nominal  motion,  and  do  not 
ignore  the  interaction  of  angular  rates  and  deflections.  They  do  assume  small  deflections  of  the  links  which 
can  be  described  by  a  summation  of  the  modal  shapes  and  a  linear  model  of  elasticity.  Only  rotational  joints 
arc  allowed.  The  results  arc  quite  tractable  for  automated  computer  solution  of  arbitrary  rotary  joints. 
Preliminary  programs  written  to  evaluate  computational  efficiency  show  that  this  method  requires  about  2.7 
times  as  many  computations  as  the  most  efficient  rigid  formulations  with  the  same  number  of  degrees  of 
freedom.  1110  rigid  model  could  incorporate  21  degrees  of  freedom  compared  to  12  degrees  of  freedom  (6  of 
which  arc  joints)  for  this  flexible  model.  Thus.  IS  degrees  of  freedom  in  the  rigid  model  could  be  used  to 
approximate  the  flexibility  that  the  6  flexible  degrees  of  freedom  of  lltc  model  presented  here  approximate. 
Ihc  relative  accuracy  of  the  two  approximations  has  not  been  determined.  These  issues  arc  discussed  in  more 
detail  in  the  Conclusions. 

2.  Flexible  Arm  Kinematics 

Ihc  previous  works  on  rigid  arm  dynamics  use  the  serial  nature  of  manipulator  arms  which  results  in 
multiplicative  terms  in  the  kinematics.  'Ihc  modal  representation  of  flexible  structure  dynamics,  on  the  other 
hand,  is  a  parallel  or  additive  representation  of  the  system  behavior.  One  of  the  contributions  of  this  paper  is 
to  resolve  this  difference  in  a  concise  way.  As  with  many  of  the  previous  works  on  rigid  dynamics,  the  4x4 
matrices  of  Denavit  and  Hartenberg  [1 1]  arc  used.  Sunada  and  Dubowsky  [31]  used  this  representation  for 
their  flexible  arm  simulations  but  did  not  produce  a  complete  nonlinear  dynamic  simulation.  Other  workers 
such  as  Hughes  [17]  relied  on  the  more  general  formulation  provided  by  a  vector-dyadic  representation. 
While  Silver  [27],  Hollcrbach  [15],  and  others  have  pointed  out  the  relative  inefficiency  of  the  4x4 
formulation,  the  conceptual  framework  is  most  advantageous  when  tackling  the  complexity  of  the  flexible 
dynamics. 

Define  the  position  of  a  point  in  Cartesian  coordinates  by  an  augmented  vector: 

[1  x-componcnt  y-componcnt  z-componcnt]T. 

Define  the  coordinate  system  [x  y  z]j  on  link  i  with  origin  0{  at  the  proximal  end  (nearest  the  base)  oriented  so 
that  the  x  axis  is  coincident  with  the  neutral  axis  of  the  beam  in  its  undeformed  condition.  Ihc  orientation  of 
the  remaining  axes  will  be  done  so  as  to  allow  effleient  description  of  the  joint  motion.  A  point  on  the  neutral 
axis  at  x  =  ij  when  the  beam  is  undeformed  is  located  at  'hj(tj)  under  a  general  condition  of  deformation,  in 
terms  of  system  i. 

By  a  homogeneous  transformation  of  coordinates  the  position  of  a  point  can  be  described  in  any  other 
coordinate  system  j  if  the  transformation  matrix  jW.  is  known.  Ihc  form  of  this  matrix  is 


.  ****** 
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where 


1  |  0T 
x.  component  of  Oj  | 
j/j  component  of  Oj  1 
Zj  component  of  O.  | 


JR  =  a  3x3  matrix  of  direction  cosines 
0  =  a  1x3  vector  of  zeros. 


Thus  in  terms  of  the  fixed  inertial  coordinates  of  the  base  the  position  of  a  point  on  link  i  is  given  as 


(1) 


hi  =  0W.ihi  =  W,^  (2) 

where  the  special  ease  of  °\Vj  =  W..  It  is  useful  to  separate  the  transformations  due  to  die  joint  from  the 
transformation  due  to  the  flexible  link  as  follows 


Wt  =  Wr.  "V.  Ai  =  Wj-iAi 


where 

Aj  =  the  joint  transformation  matrix  for  joint  j 

1^  !  =  die  link  transformation  matrix  for  link  j- 1  between  joints  j^l  andj 

=  the  cumulative  transformation  from  base  coordinates  to  O.^  at  the  distal  end  of  link  j. 


(3) 


0H  is  fixed  to  the  link  j-1  and  with  no  deflection  ( x  y  z]^,  is  parallel  to  lx  y  z]^  with  x coincident  with 
XH‘ 


To  incorporate  the  deflection  of  the  link,  the  approach  of  modal  analysis  is  used  which  is  valid  for  small 
deflection  of  the  link. 


\<1»)  = 

1 
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0 

m. 

+  E«ij 
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P  1 
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V’> 

Yi/n) 

.  0  . 

where 

Xy .  Vy,  z...  s  the  Xj.  yt,  and  z.  displacement  components  of  mode  j  of  link  fs  deflection,  respectively. 

Sj.  =  the  time  varying  amplitude  of  mode  j  of  link  i 

m,  =  the  number  of  modes  used  to  describe  the  deflection  of  link  i. 

The  link  transformation  matrix  must  also  incorporate  the  deflection  of  the  link.  Here  the  rotations  as  well  as 
the  translations  of  the  deflection  must  be  represented.  If  one  consistently  requires  small  rotations  the 
direction  cosine  matrix  simplifies  as  noted  in  (9]  and  furthermore  the  small  angles  can  be  assumed  to  add 
vcctoraliy.  'Iliis  is  basic  to  the  approach  used  here.  The  link  transformation  matrix  can  then  be  written  as 

E  =  [H.  +  E  (5) 
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Mij  = 


and  where 

AH  variables  in  brackets  arc  evaluated  at  1 

8  0  ...  8  ..  =  the  x..  y  .,  and  /..  rotation  components  of  link  i,  respectively. 

I.  =  the  length  of  link  i. 

To  find  the  velocity  of  a  point  on  link  i,  take  the  time  derivative  of  the  position: 


(7) 


dht  =  h.  =  W.’hj  +  W.'hj. 
dt  ' 


(8) 


Due  to  the  serial  nature  of  the  kinematic  chain,  it  is  computationally  efficient  to  relate  the  position  of  a  point 
and  its  derivatives  to  proceeding  members  in  the  chain.  By  differentiating  2  one  obtains: 


where 


wi“wHAJ+wi-lAi 

(9) 

00) 

_>• 

II 

. 

01) 

*i  =  tJW+uA 

(12) 

Uj  =  3A/9qj 

u2j  =  a2yaq,2 

q.  =  the  joint  variable  of  joint  j. 

Thus  Wj  and  W.  can  be  computed  recursively  from  W.  |t  its  derivatives,  and  the  partials  with  respect  to 
the  variables  of  link  j-1  and  joint  j.  No  mixed  partials  arc  explicitly  present  This  computational  approach  is 
similar  to  that  proposed  by  Hollcrbach  (1 5]  for  rigid  link  arms.  Here  one  additionally  needs  W._  (  and  its 
derivatives.  These  can  be  computed  recursively  from  W..t  and  its  derivatives: 


Wj  =  WjEj 


03) 
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Wj  =  Wj  Kj  +  Wj  Kj 

Wj  =  Wjt'  +2Wi'  +  Wj  kj 
mj 

►)  -  2  «,  M„ 

k=  l 


mj 

fi,  =  E  «jk  MJk 

k  =  I 


(14) 

(15) 

(16) 

(17) 


The  last  two  equations  illustrate  how  the  deflection  transformations  enter  even  more  simply  into  the 
kinematics  on  a  per  variable  basis  than  do  the  joint  variables.  Ihis  is  due  to  the  small  deflection  assumption 
and  the  form  chosen  for  the  transformation.  iTic  recursive  nature  of  the  velocity  and  acceleration  is  preserved 
from  the  rigid  ease.  For  the  simulation  equations  the  terms  involving  second  derivatives  of  die  joint  and 
deflection  variables  will  be  separated  from  the  above  expressions  and  included  in  die  inertia  matrix  to  make 
up  the  coefficient  matrix  of  the  derivatives  of  the  state  variables.  ITic  "inverse  dynamics"  solution  that 
proceeds  directly  from  the  lagrangc  formulation  has  little  obvious  utility. 

3.  System  Kinetic  Energy 

In  this  section  the  expression  for  the  system  kinetic  energy  is  developed  for  use  in  lagrangc’s  equations. 
First,  the  kinetic  energy  for  a  differential  element  is  written.  ITicn.  integration  of  this  differential  kinetic 
energy  over  the  link  gives  the  link’s  total  contribution.  'ITiis  produces  terms  that  arc  the  equivalent  of  the 
moment  of  inertia  matrices  of  rigid  link  arms.  Summation  over  all  the  links  provides  the  total  kinetic  energy. 


The  kinetic  energy  of  a  point  on  the  i-th  link  is 

dkj  =  _Ldm  Tr{  hjhT} 

2 


(18) 


where 


dm  is  the  differential  mass  of  the  point  and 
Tr{.}  is  the  trace  operator. 


Hxpanding  18  and  using  the  fact  that  Tr{A  Br}  =  TrjB  AT}>  the  expression  fordkj  becomes 
dk  =  _Ldm  Tr{  W.  *h,  Vr  wj  +  2  W  V  V1  +  W.  V  HtJ  W J  } 


(19) 


where 


m, 


x..  y..  /.. 
|J  I 


3T 


(20) 


V  t  W 

j=i 

By  integrating  over  the  link  one  can  obtain  the  tout  link  kinetic  energy.  In  this  report  it  is  assumed  that  the 
links  arc  slender  beams  because  it  makes  the  central  development  clearer.  Other  mass  distributions  could  be 
used  with  a  slight  departure  here  in  the  development.  For  slender  beams  dm  =  p  dij  and  one  can  integrate 
over  ij  from  0  to  1..  Only  the  terms  in  'h.  and  its  derivatives  arc  functions  of  ij  for  this  link.  Thus  the 
integration  can  be  performed  without  knowledge  of  W.  and  its  derivative.  Summing  over  all  n  links  one  finds 
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the  system  kinetic  energy  to  be 


n 

r 

K  = 

E 

I  “ki 

i=  1 

0 

n 

K  = 

E 

1  r{  W.  nywU  2Wj  B^  wT  +  Wj  Bj |  wT  } 

i=l 

where 

1. 

»  I 

B,i  = 

J 

M  \  'h'  dif. 

2 

0 

(21) 


(22) 


(23) 


lly  interchanging  the  integration  in  23  and  the  summations  involved  in  the  definition  of  V  in  20  one  obtains 

A  A  •  • 

»«i=  E  E  <24> 

j=l  k=l 

where 


i. 


-.i,  *  x / 


C-  =  -Ly  M  [0  *ik  yik  1  (0  Vij  z,jl  dn- 

2  o 


(25) 


Cjkj  has  units  of  an  inertia  matrix  and  serves  a  similar  function.  While  shown  here  as  a  4x4  matrix  it  is 
nonzero  only  in  the  3x3  (lower  right).  It  can  also  be  shown  that  Cjkj  =  C\’k.  lly  choosing  the  assumed  mode 
shapes  in  an  appropriate  manner,  it  is  possible  to  reduce  the  number  of  nonzero  terms  in  24.  This  matter  is 


discussed  in  light  of  computational  speed  in  die  conclusions. 
The  other  tejms  in  equation  22  can  similarly  be  found: 


«*  = 

J 

2  0 

fi  ‘hj'iiT  dr) 

m. 

nr  m. 

«2i  = 

£ 

*yS  +  £ 

j=l 

k  =  l  j=l 

where 

,'i 

s= 

-L / 

2  0 

[1  T|  0  0]1  [0  xy  yyZyldij. 

(26) 


(27) 


(28) 


Finally,  by  a  similar  approach: 

B  *  =  _!_/  fiVhjrdi, 


2  o 


m.  «»ij 

BJi  =  ci+  E  Vc*  +  c5  +  E  E  Myc, 


m.  m. 


j=l 


k=l  j=l 


tt 


(29) 


where 
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C ,-xJ  [1  T|00)r[l  1)00]  dij.  (30) 

2  o 

Ibis  final  term  contains  the  rigid  body  inertia  terms. 

It  should  be  noted  that  these  terms  are  easily  simplified  if  one  link  in  the  system  is  to  be  considered 
rigid,  in  which  nrv  =  0.  Should  a  link  consist  of  a  flexible  member  with  rigid  appendages  the  above  derivation 
is  readily  extended  to  modify  the  matrices  Cjkj.  Cjk.  and  C.  with  no  further  modifications  to  the  succeeding 
development.  In  fact,  these  matrices  could  be  obtained  by  finite  element  analysis  should  die  link  shape  be 
irregular  as  is  often  the  case.  Furthermore,  the  expression  for  BJ4  contains  a  term  of  order  i2  which  is  by 
definition  small  and  a  candidate  for  later  elimination.  Finally,  much  of  the  complexity  of  the  integration  of 
the  modal  shape  products  can  be  done  offline,  once,  for  a  given  link  structure. 


3.1 .  Derivatives  of  Kinetic  Energy 

For  construction  of  I  .agrange's  equations  one  needs 

3K/3a,  3K/36.f.  -^  (  3K  /  3a  )  .  and-^  (  3K  /  35.f ) 
j  dt  J  dt  J 

First  consider  3K  /  3q..  This  will  involve  the  purtials  of  all  the  terms  in  22.  some  of  which  arc  zero.  In 
fact  only  W.  for  j  <  i  <;  n  provides  nonzero  partiais  with  respect  to  q..  The  time  derivative  of  the  partial  is 
then  taken.  In  this  respect  the  following  equivalences  should  be  noted: 


3W;  /  3^  =  3WS  /  3qj  (31) 

—  (  3W  /  3q  )  =  3W  /  3q  (32) 

dt  1  J  J 

3^/38.,=  31^/35^  (33) 

j  (  3W.  /  38jf )  =  3W.  /  38jf .  (34) 


Also  helpful  in  simplifying  the  result  is  that  Tr{A}  =  Tr{AT}  for  any  square  matrix  A  and  that  B3j  is 
symmetric.  Considerable  cancellation  and  combination  results  when  the  terms  in  I  .agrange's  equation 
involving  the  kinetic  energy  are  combined.  The  result  of  this  combination  is 

|(3K/3qj)-3K/3qj  = 


3  W; 


^  e  t  «»<cil+ci+ 1  «)]*; 


l=J 

m. 


k  =  1 


1=1 


+  C  £  li  «iCB»J'+[2  2;  «,(C,+  £  «acJ]wT]} 

1=1 


k  =  l 


k  =  1 


1=1 


(35) 


I 

s' 


g 


Note  above  terms  of  the  form  6jk  6j(  which  arc  second  order.  These  can  be  ignored  consistent  with  the 
assumption  that  the  deflections  arc  small.  Noting  the  rccurrancc  of  certain  terms  above,  it  is  convenient  to 
define  the  following: 


-  c,.  +  E 

1  =  1 


k=  1 


(36) 

(37) 


When  these  definitions  arc  substituted  into  equation  35  one  obtains: 


-(  dK/dqj-d  K/9q.  = 
dt  ‘  ‘ 


n  aw  "i  ™L 

2  E  E  K  >>*  w’+2  E 

i=j  oq:  k=l  k  =  1 


(38) 


The  partials  of  K,  with  respect  to  5.f  and  arc  considerably  more  complex  due  to  the  fact  that  R)jt  B2j, 
and  B3j  arc  functions  of  the  deflection  variables.  The  techniques  of  simplification  arc  similar.  An  additional 
simplification  arises  due  to  the  fact  that  if  A  were  any  antisymmetric  matrix,  and  if  W  were  a  matrix 
compatible  for  multiplication,  then  Tr{  W  A  W1  }  =  0.  An  antisymmetric  matrix  occurs  from  the  difference 
of  a  matrix  and  its  transpose. 

^(3K/3fijr)-3K/a«jr  = 


2  £  Tr(4^ 


nr  nr 

[G,W ?  +  £  ^DbWT  +  2  £  + 


i=j  +  l  “'’it  k=  1  k  =  1 

Tr{2[  W,VJ*I  2-  *»Sr+Wi  2- 


k  =  l 


k=l 


(39) 


4.  System  Potential  Energy 

The  potential  energy  of  the  system  arises  from  two  sources:  clastic  deformation  and  gravity.  In  both 
eases  they  arc  included  by  first  writing  die  potential  energy  contribution  of  a  differential  element,  integrating 
over  the  length  of  the  link,  and  then  summing  over  all  links. 


4.1 .  Elastic  Potential  Energy 

Consider  a  point  on  the  i-th  link  undergoing  small  deflections.  First  restrict  the  link  of  the  slender 
beam  type.  The  clastic  potential  is  accounted  for  to  a  good  approximation  by  bending  about  the  transverse  y, 
and  /.)  axes  and  twisting  about  the  longitudinal  x.  axis.  Compression  is  not  initially  included  since  it  is 
generally  much  smaller.  Along  an  incremental  length  dtj  the  clastic  potential  is 
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$xi,  Oyi,  and  $fi  arc  the  rotations  of  the  neutral  axis  of  the  beam  at  the  point  ij  in  the  xr  and  z, 
directions,  respectively.  Since  deflections  arc  small,  these  directions  arc  essentially  parallel  or 
perpendicular  to  the  neutral  axis  of  die  beam. 

K  =  Young's  modulus  of  elasticity  of  the  material 

G  =  Ihc  shear  modulus  of  the  material 

lx  =  Ihc  polar  area  moment  of  inertia  of  the  link  cross  section  about  die  neutral  axis. 

Iy,  1^  =  the  area  moment  of  inertia  of  the  link  cross  section  about  the  yt  and  r  axes,  respectively. 

With  a  truncated  modal  approximation  for  the  beam  deformation  the  angles  <?xj,  0yj,  and  0 A  are 
represented  as  summations  of  modal  coefficients  times  the  deflection  variables.  Ihc  x  rotation,  for  example  is 


8  =  >  ..  . 

xi  £—1  ikxik 


where  0xA  is  the  angle  about  the  xi  axis  corresponding  to  the  k-th  mode  of  link  i  at  the  point  n.  When  dvd  is 
integrated  over  the  link  the  integration  can  be  taken  inside  the  modal  summations  of  equation  41  and  its 
corresponding  y  and  /.  components.  Ihc  following  definitions  then  prove  useful: 

Kikl  =  Kxikl  +  Kyikl  +  Krikl  •  <42> 


where 


l. 

xikl  =  l 


G 


■ 


d$  ,  de 


-t 


do dd  _ 


Note  that  Kjk|  =  K.|k  and  that  for  certain  special  eases  the  orthorgonality  of  the  modal  functions  can 
eliminate  many  of  the  terms  in  equations  43, 44,  and  45.  Ihc  clastic  potential  for  the  total  system,  Vc  can  then 
be  written  as 


m.  m. 


■■  ii«<  in? 

=  -J-  £  E  E  Ww 

2  i=i  k  =  l  1=1 


(46) 


(47) 


Note  that  the  V#  is  independent  of  the  joint  variables. 


3V, 

3qi 


=  0. 


Kor  deflection  variables 

av  A 

ds  r  ~  sik  Sr  (48) 

jf  k=  1 

llic  form  of  equation  48  is  much  more  general  than  the  initial  assumptions  made  regarding  the  contributions 
to  the  clastic  potential  energy  would  allow.  Compression  strain  energy,  and  link  forms  other  than  beams  can 
be  represented  in  this  form.  lire  values  of  the  coefficients  K.^can  be  determined  analytically  or  numerically, 
eg.  by  finite  element  methods. 


4.2.  Gravity  Potential  Energy 

For  a  differential  element  on  the  i-th  link  of  length  dij  the  gravity  potential  is 

dvgi  =  -/i8rW.ihjdi,. 

where  the  gravity  vector  g  has  the  form 
*  ‘  =  1 0  g»  gy  g*  1  • 

When  integrated  over  the  length  of  the  beam  and  summed  over  all  beams,  the  gravity  potential  becomes 


(49) 


V.=  «T  E  wi', 

i= 1 


where 


+  £  *»•» 

k=i 

A/.  =  the  total  mass  of  link  i 

rri  =  1 1  r„  0  0  J,  a  vector  to  the  center  of  gravity 
from  joint  i  (undeformed) 

«ik  =  /  f*I0  xjk  ytt  z.k]Tdi|. 


(50) 


(51) 


(52) 


Note  that  ejk  is  found  in  the  top  row  of  Cjk.  It  is  the  distance  from  the  undeformed  center  of  gravity  to  the 
center  of  gravity  when  all  6  arc  zero  except  5jk .  which  is  one.  Ihc  total  distance  to  the  center  of  gravity  from 
O.  (joint  i)  is  multiplied  by  the  mass  to  give  r(. 

Upon  taking  the  partial  derivatives  required  by  Lagrange’s  equations  we  find  fbr  the  joint  variables 


T 


u 


ave  A  aw. 

For  the  deflection  variables,  for  1  £  j  $  n-1 
3V  "  /  3W;  . 

34if  i=j  +  l  34jf 


«1wjV 


Forj  =  n 
3V, 
36 


=  •8' Wne 


nf 


n  nf 


(53) 


(54) 


(55) 


5.  Lagrange’s  Equations  in  Simulation  Form 

At  this  juncture  the  components  of  the  complete  equations  of  motion  in  I  ^grange's  formulation,  except 
for  the  external  forcing  terms,  have  been  evaluated  in  equations  38.  47.  and  S3  for  the  joint  equations;  and  in 
equations  39.  48.  54  and  55  for  deflection  equations.  The  external  forcing  terms  arc  the  generalized  forces 
corresponding  to  the  generalized  coordinates:  the  joint  and  deflection  variables  in  this  ease.  The  generalized 
force  corresponding  to  joint  variable  q,  is  the  joint  torque  F.  For  the  deflection  variables  the  corresponding 
generalized  force  will  be  zero  if  the  corresponding  modal  deflections  or  rotations  have  no  displacement  at 
those  locations  where  external  forces  arc  applied.  Thus  it  is  assumed  for  the  present  development  that  the 
modal  functions  arc  selected  so  that  is  the  ease.  This  is  convenient  for  using  the  results  as  well.  All  motion  at 
the  joint  is  described  in  terms  of  the  joint  variable.  (This  is  not  true  in  the  approach  taken  by  Sunada  and 
Dubowski  [31].)  The  form  of  1  ^grange's  equations  will  then  be; 


The  joint  ttwatiunj 

d  ,  .  .  3V  3V 

-  (3K/3q.)  -  3K/3q.  +  — -*■  +  — =  F,  . 
dt  3q  3q  i 


(56) 


TO£-d£fltttiflqgq»atifln  j.f 

d  .  .  .  3V  3V 

-  (3K/3«jr )  •  3K/Mjr  +  *  0  ■ 


(57) 


These  equations  are  in  the  "inverse  dynamic"  form.  I'o  convert  them  to  the  simulation  form  one  must  extract 
the  coefficients  of  the  second  derivatives  of  the  generalized  coordinates  to  compose  an  inertia  matrix  for  the 
system.  The  second  and  first  derivatives  together  make  up  the  derivative  of  the  state  vector,  which  can  be 
used  in  one  of  the  available  integration  schemes,  c.g.  Runga-Kutta.  to  solve  for  the  state  as  a  function  of  time 
forgiven  initial  conditions  and  inputs  F(. 


5.1.  Kinematics  Revisited 

Hie  purpose  of  this  section  will  be  to  extend  the  kinematics  to  separate  the  second  derivatives  of  the 
joint  variables  and  deflection  variables  from  the  expressions  for  W.  and  W..  Other  occurrences  of  these 
derivatives  arc  already  explicit  in  the  formulation  as  it  exists. 


First  consider  the  product  of  transformations  which  make  up  W.  and  two  alternative  ways  of  expressing 


it 
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Wj  =  MiV^-W-Mi 


Carrying  through  the  derivatives  one  obtains 


i 


E  (  *IH  Uh  "  Wi‘>h  + 

h  =  1 


E  ".“v"*!  *«>  +  *.. 

k  =  1 


(58) 

(59) 


(60) 


For  the  corresponding  expression  for  W.  write 


=  \A“*. 

-  w»  *wi  • 


EAiU.*** +  E  2»i 


(61) 

(62) 

(63) 


h=l  h=l  k= 1 

X  „  X 

The  value  of  Wy.  and  Wyj  can  be  calculated  recursively  as  shown  in  equations  15  and  10.  respectively,  for  VV 
and  W.  by  only  eliminating  terms  involving  (jj  and  5jk.  The  result  is 

+  +  (64) 


Wyj  =  WyjEj  +  2WjEj. 


(65) 


5.2.  Inertia  Coefficients 

To  obtain  the  inertia  coefficients  that  multiply  the  second  derivatives,  substitute  equations  63  and 
60  into  the  relevant  parts  of  the  equations  of  motion,  equations  38  and  39.  respectively.  Collecting  the  terms 
and  arranging  them  for  efficient  computation  requires  the  steps  outlined  in  this  section. 

5.2.1 .  Inortia  Coefficiants  of  Joint  Variables  in  the  Joint  Equations 

All  occurnnccsofq.  in  equation  38  arc  in  the  expression  for  W.1.  When  these  terms  arc  isolated,  a 
double  summation  over  tnc  indices  i  and  h  exists.  Interchange  the  order  of  the  summation  as  follows: 
n  i  n  n 

E  E  *  E  E  • 

i=j  h  =  1  h  =  l  i=max(h.j) 

The  resulting  coefficient  for  joint  variable  qh  in  the  joint  equation  j  is 
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where 

n 

jph=  E  j^iGih^ir  (67) 

i=max(h.j) 

Note  that  if  one  exchanges  j  and  h  and  transposes  inside  the  trace  operation  an  identical  expression  is 
obtained.  Iliis  indicates  the  symmetry  of  the  inertia  matrix  which  is  used  to  reduce  the  number  of 
computations  required.  1110  expression  foHKh  can  be  computed  recursively:  this  will  be  described  later  to 
further  improve  the  efficiency  of  calculation. 


5.2.2.  Inertia  Coefficients  of  the  Deflection  Variables  in  the  Joint  Equations 

ITic  deflection  variables  appear  both  in  the  expression  for  W 1  and  explicitly  in  equation  38.  After 
substituting  W('  into  equation  38.  collect  terms  in  fi.r  and  exchange  Lite  order  of  summations  as  follows 
n  i-1  n-1  n 

E  E  -  E  E 

i=j  h  =  1  h  =  l  i  =  max(h  +  l,j) 

The  resulting  coefficient  of  5hk  in  joint  equation  j  is  Jjhk-  The  terms  to  be  included  depend  on  the  relative 
values  of  j  and  h.  Ilic  following  hold  for  1  ^  k  ^  mh> 

Forh  =  n,  j  =  1 ...  n: 


(68) 

forh  =  j...  n-l.j  =  1  ...n-1: 

V  -2Tr{(WhlUj)[  lK„  »£+  <*„  IV  ]  wj  >; 

(69) 

forh  =  1 ...  j-1,  j  =  2...  n: 

(70) 

where  for  h  =  1 ...  n-1,  j  =  1 ...  n: 

V  E 

i=max(h  +  l,j) 

(71) 

It  can  be  shown  that  the  inertia  coefficient  for  the  deflection  variable  {>k  in  the  joint  equation  j  is  the 
same  as  the  coefficient  for  the  joint  variable  q^  in  the  deflection  equation  h.k.  This  further  extends  the 
symmetry  of  the  inertia  matrix  and  reduces  the  necessary  computation; 


5.2.3.  Inertia  Coefficients  of  the  Deflection  Variables  in  the  Deflection  Equation 

In  a  manner  similar  to  the  previous  two  types  of  coefficients,  the  inertia  coefficients  of  the  deflection 
variables  in  the  deflection  equations  arc  evaluated.  Symmetry  of  the  coefficients  can  be  shown  such  that  the 
coefficient  of  variable  h.k  in  equation  j.f  is  the  same  as  the  coefficient  of  variable  j,f  in  equation  h,k. 
Substituting  equation  63  into  equation  39,  isolating  the  second  derivatives  of  the  deflection  variables,  and 
interchanging  the  order  of  summations  enables  the  inertia  coefficients  to  be  identified.  Further  simplification 
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is  based  on  the  identity  that,  for  any  three  square  matricics  A,  B  and  C 
Tr{A  BC}  =  TrfCA  B}  =  Tr{BC  A}. 

Furthermore  die  rotation  matrices  in  the  transformation  matrices  arc  orthogonal  so  that  R(  R  ®  =  I,  a  3x3 
identity  mauix.  This  coupled  with  die  zero  first  row  and  column  of  C.kf  results  in  an  especially  simple  form 
for  two  of  the  four  eases.  I  he  following  hold  for  1  <  k  ^  mh  and  1  £  f  s  ny 
Forj  =  h  =  n: 

Forj  =  h  =1 ...  n-1: 

W  =  rlr{MAMi+cj.f>-  ,73) 

Forh  =  n;  j  =  1 ...  n-1: 

W  =  2l'{w,M,r <«> 

Forj  =  1 ...  n-1:  h  =  j+1 ...  n-1: 

■ju  -  2  1'H  [V  M7  +  iw „  ■>„]«£}  . 

(75) 

Terms  in  the  above  defined  forj  =  1 ...  n-1 :  h  =  1 ...  n-1  arc: 
n 

j*h=  £  ^  G,  hW®  .  (76) 

i=max(j+l.h  +  l) 

5.2.4.  Recursions  in  the  Calculation  of  the  Inertia  Coefficients 

Since  the  inertia  matrix  is  a  square  matrix  it  requires  the  calculation  of  n*  terms  where  n(  is  the  total 
number  of  variables: 

n 

n,  =  n  +  £  mi 
i=  1 

The  fact  that  the  matrix  is  symmetrical  reduces  the  number  of  distinct  terms  to  n((nt+ 1)/2.  which  still  has  a 
second  power  dependence.  Thus  while  the  inverse  dynamics  computation  complexity  can  be  made  linear  in 
n(,  simulation  requires  the  inertia  matrix  with  complexity  dependent  on  n*.  Since  n,  can  be  quite  large  for 
practical  arms  it  is  important  to  reduce  the  coefficient  of  the  squared  term  as  much  as  possible.  I)uc  to  their 
short  or  even  /.cro  length,  it  is  possible  for  some  links  to  be  essentially  rigid.  Anthropomorphic  arms,  for 
example,  have  two  links  which  arc  much  longer  than  the  others  and  tend  to  dominate  the  compliance.  Many 
of  the  terms  derived  above  may  not  be  needed  for  these  links,  four  of  the  six  links  in  the  anthropomorphic 
example.  Any  recursive  scheme  fur  calculating  the  terms  in  the  equations  should  not  require  these 
calculations  as  a  means  to  get  to  needed  terms. 

Consider  the  calculation  of  equations  67,  71,  and  76.  Several  recursive  schemes  could  be  arranged  for 
the  efficient  calculation  of  these  quantities.  Fquation  71  is  only  needed  if  the  link  corresponding  to  the 
variable,  link  h,  is  flexible.  That  is,  if  mh  >  0.  liquation  76  is  only  needed  if  both  the  link  of  the  variable  and 
the  link  of  the  equation,  link  j,  is  also  flexible.  Thus  we  propose  the  following  recursive  scheme  for 


T- 

✓ 


♦ 


IS 


calculating^,  jFh.  and M>h.  'rhc  following  hold  for  1  £  k  s  mh;  I  £  f  £  ntj . 
Initialization: 


6 

ll 

(77) 

Kor  j  >  h  s  n: 

(78) 

Kor  j  =  h: 

-  ch  + 

(79) 

lfmh  >  0 calculate: 

it-  =  if  at 
rh  rh+l/V 

(80) 

If  mh  >  0  and  m.  >  0  calculate: 

laf 

+ 

<r 

ii 

(81) 

5.3.  Assembly  of  Final  Simulation  Equations 

The  complete  simulation  equations  have  now  been  derived.  It  remains  to  assemble  them  in  final  form 
and  to  point  out  some  remaining  recursion  relations  that  can  be  used  to  reduce  the  number  of  calculations. 
The  second  derivatives  of  the  joint  and  deflection  variables  arc  desired  on  the  "led  hand  side”  of  the  equation 
as  unknowns  and  the  remaining  dynamic  effects  and  the  inputs  arc  desired  on  die  "right  hand  side."  To  carry 
out  this  process  completely  one  would  take  the  inverse  of  the  inertia  matrix  J  and  prcmultiply  the  vector  of 
other  dynamic  effects.  This  inverse  can  only  be  evaluated  numerically  because  of  its  complexity.  Ill  us  for  the 
present  purposes  the  equations  will  be  considered  complete  in  the  following  form: 

J  i  -  R,  (82) 


where 


J  =  Inertia  matrix  consisting  of  coefficients  previously  defined  in  the  order  for  multiplication 
appropriate  for  z 

z  =  the  vector  of  generalized  coordinates 
=  1<1|  *12  -  *|mj  ^2  *21—  *2m2 “ qh  *hl  ~  ~  ~  ®nmj 

qh  =  the  joint  variable  of  the  h-th  joint 

*hk  =  the  deflection  variable  (amplitude)  of  the  k-thmodc  of  link  h 


R  s  vector  of  remaining  dynamics  and  external  forcing  terms 
*  fRl  RU  RI2  “•  Rlm1  R2  R21  ~  R2m2  "  Rj  RJI  “  Rjf  ~  Rjm^ 
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Rj  =  dynamics  from  the  joint  equation  j  (equation  S6)  excluding  second  derivatives  of  the 
generalized  coordinates 

Rjr  =  dynamics  from  the  deflection  equation  jf  (equation  57)  excluding  second  derivatives  of  the 
generalized  coordinates 


The  elements  of  J  have  just  been  formulated  and  can  be  arranged  to  form  the  proper  equations  in  the  order 
described  above.  Ibis  order  has  been  selected  because  it  results  in  the  symmetrical  appearance  of  J.  The 
elements  of  K  have  not  been  explicitly  given  with  the  second  derivatives  removed.  Ihcsc  arc  given  below 
with  some  recursions  to  facilitate  their  compulation. 

R,  =  *2Tr  {  U,  Q,  }+  gTU,  P,  +  F,  (83) 


R  =  -2lr  {  U  Qj  >  +iT*J.1UJ«i+FJ 

mn  tn 

E**c.»rlwI>-  E  K„,+ s' 'V„f 


k=l 


k  =  1 


k  =  l 


RJf  =  -2  Tr  {  Wj Mj,  V ,  Q,. ,[  W,j  l)Jf  +  2 W,  £  «, Cjl( ]  W'  > 

•  E  ®jk Kjir+  *TwiMjrAj+i,*j+i+  *Twi*jr 


k=l 


where 


m. 


Q.  «  C.  K*  2  (  E  «» ,■>*)*! 

k  =  l 

Qj  =  Gj wI + 2 (  £  *jkVw/+EiViQi+i 

k=l 

=  E  «„>  •„> 

k=l  • 

A 

P.  =  A/,  r.  +  ^2  *jk  +  Aj+ 1  Pj+ 1 
k  =  l 


(84) 

(85) 


(86) 

(87) 

(88) 

(89) 

(90) 
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6.  Conclusions 

The  above  model  is  successful  in  terms  of  its  accuracy  and  its  speed,  lire  two  qualities  arc  somewhat 
related  in  that  accuracy  of  the  flexible  representation  can  be  improved  by  increasing  the  number  of  modes 
used  to  represent  the  link  deflection  at  the  expense  of  calculation  time.  'ITic  issue  is  further  complicated  by 
the  choice  of  mode  shapes,  range  of  motion  considered,  and  the  arm  configuration.  Furthermore,  limited 
information  is  available  in  the  literature  for  comparison.  A  simple  comparison  has  been  used  in  the  past  and 
can  be  performed  for  calculation  complexity.  Hollcrbach  (IS]  compares  several  approaches  to  the  inverse 
dynamics  problem  of  rigid  arms  by  different  authors.  Walker  (33(  gives  a  similar  count  for  four  approaches  to 
the  simulation  problem.  Sunada  (31]  has  given  computation  times  for  a  given  manipulator,  trajectory,  and 
computer  for  his  flexible  simulation.  Comparison  to  the  calculation  counts  of  rigid  models  arc  given  for  a 
rough  comparison  of  speeds  in  this  section.  No  attempt  at  a  quantitative  comparison  of  the  accuracy  is  made. 

To  determine  the  number  of  calculations  from  the  equations,  a  choice  must  be  made  on  how  some 
matrix  products  arc  implemented.  Hollcrbach  chose  to  use  the  most  straightforward  implementation  of  the 
equations.  The  approach  here  is  quite  different.  Obvious  simplifications  in  the  multiplication  of  matrices 
with  known  constant  rows,  the  top  row  of  a  transformation  matrix  for  example,  arc  assumed  in  these 
computations.  Hie  4x4  matrix  transformation  was  chosen  for  its  conceptual  convenience  and  die  calculation 
count  will  not  be  intentionally  penalized  for  that  choice.  Furthermore,  certain  products  appear  in  multiple 
equations  and  arc  assumed  to  be  saved  when  needed  later.  Special  purpose  multiply  routines  arc  used 
whenever  they  can  capitalize  on  die  special  structure  of  a  given  matrix.  Finally,  in  the  simulation  form  the 
calculations  needed  to  invert  the  inertia  matrix  arc  not  included,  and  no  consideration  is  given  to  the 
calculations  of  the  integration  routine.  The  general  form  of  the  modal  parameters  arc  used  however.  This 
results  in  all  combinations  of  modes  h  and  k  in  the  matrix  Cjhk  to  be  computed  and  used  and  hence  introduces 
a  squared  dependence  on  the  number  of  modes  on  each  inertial  coefficient  of  the  deflection  variables.  With 
these  assumptions  the  number  of  calculations  is  approximate: 

bjumfacr  uf  multiplications; 

6n|m2+  17.5  nf  m2  +  118  n2  m+  74  n  nfm  + 

137.5  nfm  +  84  n2  +  86  n  nf  +  279  n+  126  nf*  57 


Number  of  additions; 


6.5  n2  m2  +  19  nfm2  +  115.5  n2m  +  68  n  nfm  + 
123  nfm  +  85  n2  +  80  n  nf  +  329  n  +  111  nf*91 

where;  n  =  total  number  of  joints 

nf  =  number  of  flexible  links 
m  =  number  of  modes  describing  each  flexible  link 


1110  above  approximation  assumes  an  "average"  joint  complexity  over  two  common  types  of  rotary  joints,  the 
same  number  of  modes  on  each  flexible  link,  a  rigid  last  link  and  a  flexible  first  link. 

If  assumed  mode  shapes  arc  restricted  so  that  the  shape  functions  in  the  x.  y.  and  z  directions  arc 
orthogonal,  only  Cjkk  will  be  non*/.cro.  'Ihis  is  a  stronger  requirement  than  the  orthogonality  of  the  set  of 
complete  mode  shapes,  but  would  often  be  realized  with  simple  mode  shapes.  It  has  not  been  determined  if 
this  would  improve  the  combination  of  speed  and  accuracy. 

Ihis  calculation  count  can  be  roughly  compared  to  rigid  link  results  available  in  the  literature 
mentioned  above.  For  a  12  degree  of  freedom  rigid  problem  the  inverse  3x3  transformation  matrix 
formulation  requires  2.66  times  as  many  multiplies  as  the  Ncwton-Kulcr  formulation.  Walker's  method  3  (his 
best)  for  simulation  requires  4.491  multiplies.  For  6  joints,  and  two  flexible  links  with  3  modes  each  the 
method  of  this  paper  requires  approximately  12,009  multiplies.  1110  ratio  of  these  simulation  methods  is  2.67. 
almost  exactly  the  same  as  for  the  inverse  dynamic  methods  with  the  same  number  of  degrees  of  freedom.  A 
modal  representation  of  flexibility  would  be  much  more  accurate  than  adding  6  imaginary  joints  to  represent 
compliance,  but  one  could  expect  to  use  15  imaginary  joints  and  6  real  joints  with  Walker's  method  with 
fewer  multiplies  than  with  the  method  of  this  paper. 

Thus  it  seems  that  in  order  to  be  competitive  with  possible  Ncwton-F.ulcr.  non-transfer  matrix 
approaches,  the  simplification  of  the  assumed  mode  shapes  will  have  to  be  made.  It  is  not  clear  that  the 
conceptual  convenience  of  the  transformation  matrix  approach  can  be  justified  relative  to  vector  dyadic 
approaches  of  Hughes  [17]  and  l.ikins[28].  Unfortunately,  computation  counts  arc  not  available  for  that 
work. 
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